Example 1: A fixed point attractor in $\mathbb{R}^2$

$F(x,y) = (\dot x, \dot y) = a(x-p_1), b(y-p_2)$



In [141]:

    
dt = 0.001

import nengo
model = nengo.Network()

# use a sigmoid neuron model
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()

# use a synapse model that only produces an output for one time step
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)

# use an exact solver that does not worry about noise
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)

with model:
    # create the neurons
    a = nengo.Ensemble(n_neurons=30, dimensions=2)
    
    # make recurrent connections that approximate this function
    def f(x):
        a = -10
        b = -10
        p1 = 0.2
        p2 = 0.3
        return dt*a*(x[0]-p1), dt*b*(x[1]-p2)    
    nengo.Connection(a, a, function=f)
    
    # add recurrent connections to turn neurons into integrators
    nengo.Connection(a, a)
    
    # record output data
    probe = nengo.Probe(a)
    
    # put in an initial stimulus to start near (-1,-1)
    nengo.Connection(nengo.Node(lambda t: [-1,-1] if t==0.001 else [0,0]), a)
    
# run the model    
sim = nengo.Simulator(model, dt=dt)
sim.run(1)

# plot the results
subplot(1,2,1)
plot(sim.trange(), sim.data[probe])
subplot(1,2,2)
plot(sim.data[probe][1:,0], sim.data[probe][1:,1])
ylim(-1,1)
xlim(-1,1)
show()

Example 2: A limit cycle in $\mathbb{R}^2$

$\dot x = -y + x(r^2 - x^2 - y^2)$

$\dot y = x + y(r^2 - x^2 - y^2)$



In [144]:

    
dt = 0.01

import nengo
model = nengo.Network()
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)
with model:
    a = nengo.Ensemble(n_neurons=30, dimensions=2)
    def f(x):
        r = 1
        return (dt*(-x[1] + x[0]*(r**2 - x[0]**2 - x[1]**2)),
                dt*( x[0] + x[1]*(r**2 - x[0]**2 - x[1]**2)))
    
    nengo.Connection(a, a, function=f)
    nengo.Connection(a, a)
    probe = nengo.Probe(a)
        
sim = nengo.Simulator(model, dt=dt)
sim.run(30)
subplot(1,2,1)
plot(sim.trange(), sim.data[probe])
subplot(1,2,2)
plot(sim.data[probe][1:,0], sim.data[probe][1:,1])
ylim(-1,1)
xlim(-1,1)
show()

Example 3: The Van der Pol oscillator

$\ddot x + \mu(x^2-1) \dot x + \omega^2 x = 0$

changing variables,

$\dot x = y$

$\dot y = -\mu(x^2-1)y - \omega^2 x$



In [146]:

    
dt = 0.1

import nengo
model = nengo.Network()
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)
with model:
    a = nengo.Ensemble(n_neurons=100, dimensions=2)
    def f(x):
        mu = 0.3
        w = 0.3
        dx = x[1]
        dy = -mu*(x[0]**2-1)*x[1] - w**2 * x[0]
        
        return (dx*dt, dy*dt)
    
    nengo.Connection(a, a, function=f)
    nengo.Connection(a, a)
    probe = nengo.Probe(a)
        
sim = nengo.Simulator(model, dt=dt)
sim.run(300)
subplot(1,2,1)
plot(sim.trange(), sim.data[probe])
subplot(1,2,2)
plot(sim.data[probe][1:,0], sim.data[probe][1:,1])
ylim(-2,2)
xlim(-2,2)
show()

Example 4: The forced Van der Pol oscillator

Basic system:

$\dot x = 0.7y + 10 x (0.1-y^2) $

$\dot y = -x$

Forcing system (whose output should be added to $\dot y$) is taken from example 2



In [184]:

    
dt = 0.01

import nengo
model = nengo.Network()
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)
with model:
    a = nengo.Ensemble(n_neurons=100, dimensions=2)
    def f(x):
        dx = 0.7*x[1]+10*x[0]*(0.1-x[1]**2)
        dy = -x[0]        
        return (dx*dt, dy*dt)
    
    nengo.Connection(a, a, function=f)
    nengo.Connection(a, a)
    probe = nengo.Probe(a)
    
    # create the forcing system
    b = nengo.Ensemble(n_neurons=100, dimensions=2)
    def f(x):
        r = 1
        rate = 1.57  # set the rate of oscillation
        return (dt*rate*(-x[1] + x[0]*(r**2 - x[0]**2 - x[1]**2)),
                dt*rate*( x[0] + x[1]*(r**2 - x[0]**2 - x[1]**2)))
    
    nengo.Connection(b, b, function=f)
    nengo.Connection(b, b)
    
    # connect the two systems together so that the first output from b is fed into the
    # second input in a
    nengo.Connection(b[0], a[1], transform=dt*0.25)
        
        
sim = nengo.Simulator(model, dt=dt)
sim.run(500)
subplot(1,2,1)
plot(sim.trange(), sim.data[probe])
subplot(1,2,2)
plot(sim.data[probe][1:,0], sim.data[probe][1:,1])
ylim(-1,1)
xlim(-1,1)
show()

Example 5: The forced Duffing oscillator

$\eta^2 \ddot x + 2 \eta \zeta \dot x + x + \alpha x^3 = f_0 cos(t)$

rearranging,

$\dot y = (-2 \eta \zeta y - x - \alpha x^3 + f_0 cos(t)) / \eta^2 $

$\dot x = y $

removing forcing term $f_0 cos(t) / \eta^2$:

$\dot y = (-2 \eta \zeta y - x - \alpha x^3) / \eta^2 $

$\dot x = y $



In [195]:

    
dt = 0.01

eta = 2.0
zeta = 0.1
alpha = 0.05
f0 = 2.5

import nengo
model = nengo.Network()
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)
with model:
    a = nengo.Ensemble(n_neurons=300, dimensions=2)
    def f(x):
        dy = (-2 * eta * zeta * x[1] - x[0] - alpha * x[0]**3) / eta**2
        dx = x[1]
        return (dx*dt, dy*dt)
    
    nengo.Connection(a, a, function=f)
    nengo.Connection(a, a)
    probe = nengo.Probe(a)
    
    # create the forcing system
    b = nengo.Ensemble(n_neurons=100, dimensions=2)
    def f(x):
        r = 1
        rate = 1  # set the rate of oscillation
        return (dt*rate*(-x[1] + x[0]*(r**2 - x[0]**2 - x[1]**2)),
                dt*rate*( x[0] + x[1]*(r**2 - x[0]**2 - x[1]**2)))
    
    nengo.Connection(b, b, function=f)
    nengo.Connection(b, b)
    
    # connect the two systems together so that the first output from b is fed into the
    # second input in a
    nengo.Connection(b[0], a[1], transform=dt*f0/eta**2)
        
        
sim = nengo.Simulator(model, dt=dt)
sim.run(100)
subplot(1,2,1)
plot(sim.trange()[2000:], sim.data[probe][2000:])
subplot(1,2,2)
plot(sim.data[probe][2000:,0], sim.data[probe][2000:,1])
ylim(-2,2)
xlim(-2,2)
show()

Example 6: Chaotic attractors in $\mathbb{R}^3 $

The Rossler system

$ \dot x = -y - z $

$ \dot y = x + a y $

$ \dot z = b + z (x-c) $



In [215]:

    
dt = 0.01

import nengo
model = nengo.Network()
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)
with model:
    a = nengo.Ensemble(n_neurons=300, dimensions=3, radius=10)
    def f(x):
        a = 0.1
        b = 0.1
        c = 9
        dx = -x[1] - x[2]
        dy = x[0] + a*x[1]
        dz = b + x[2]*(x[0]-c)
        return (dx*dt, dy*dt, dz*dt)
    
    nengo.Connection(a, a, function=f)
    nengo.Connection(a, a)
    probe = nengo.Probe(a)
        
sim = nengo.Simulator(model, dt=dt)
sim.run(300)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(sim.data[probe][:,0], sim.data[probe][:,1], sim.data[probe][:,2])
show()

The Lorentz system

$ \dot x = \sigma (y - x) $

$ \dot y = x (\rho - x) - y $

$ \dot z = x y - \beta z $



In [222]:

    
dt = 0.01

import nengo
model = nengo.Network()
model.config[nengo.Ensemble].neuron_type=nengo.Sigmoid()
model.config[nengo.Connection].synapse = nengo.Lowpass(dt)
model.config[nengo.Connection].solver=nengo.solvers.Lstsq(rcond=0)
with model:
    a = nengo.Ensemble(n_neurons=300, dimensions=3, radius=50)
    def f(x):
        sigma = 10
        beta = 8.0/3
        rho = 28
        dx = sigma*(x[1]-x[0])
        dy = x[0]*(rho-x[2])-x[1]
        dz = x[0]*x[1]-beta*x[2]
        return (dx*dt, dy*dt, dz*dt)
    
    nengo.Connection(a, a, function=f)
    nengo.Connection(a, a)
    probe = nengo.Probe(a)
        
sim = nengo.Simulator(model, dt=dt)
sim.run(50)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(sim.data[probe][:,0], sim.data[probe][:,1], sim.data[probe][:,2])
show()



In [ ]: